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INTRODUCTION 

The NASA Lewis Research Center has developed a steady state heat pipe 
code for use In heat pipe design and analysis of laboratory data, suita- 
ble for microcomputers as well as larger machines.' It has become appar- 
ent that the vapor flow algorithms in this and similar codes^-A are not 
very satisfactory at conditions where compressibility is important. Heat 
pipes having axially varying heat input, sequential evaporators, or oper- 
ating near the sonic limit, are not handled well in these codes. Typical 
of the problems In some codes is the prediction of a negative absolute 
pressure within an evaporator with large heat input, because of improper 
treatment of compressibility effects. The ability to handle compressi- 
bility and variable heat input is desired for the NASA Lewis code. There- 
fore, an accurate vapor flow algorithm concise enough to operate on a 
microcomputer at a reasonable speed and incorporating these features was 
sought. Finite difference techniques which have been used effectively 
in fundamental studies of heat pipe vapor dynamics 4 * 5 were precluded by 
code size. Among sources consulted in an attack on this problem was the 
book by Ivanovskii et al. 5 which reviews several studies leading to one- 
dimensional algorithms for compressible vapor flow in heat pipes. 

For the present work, Busse's 7 paper on pressure drop In the vapor phase 
was used as a starting point. His equations for Incompressible flow were 
readily rederived to include change in density. In addition, the usual 
assumption of a saturation temperature across each section of the pipe, 
determined by the local static pressure, was eliminated. In view of the 
studies of DeMichele, 4 a mean vapor temperature differing from the wall 
temperature was found from a boundary-layer-type energy equation. The 
resulting equations were derived so as to be suitable for liquid metal 
working fluids having two vapor species as well as for fluids with a 
single vapor species. The algorithm is Intended only for the evaporator 
and adiabatic sections of the heat pipe, because vapor flow in the con- 
denser is complicated by the onset of turbulence and reverse flow. 



ANALYSIS 


In the analysis to follow, the appropriate features of Busse's 7 method 
will be Introduced, and the development will proceed from there. In 
addition, some major assumptions for the heat pipe sections being consid- 
ered should be Indicated at the outset: Flow Is laminar. The ideal gas 

law applies to components of the vapor. The equations allow for local 
chemical equilibrium among species in the vapor. Phase equilibrium 

occurs only at the wall, so liquid droplets are absent in the vapor. f 

Following Busse, 7 no tangential velocity component exists at the wall in 
mass addition or removal, the radial pressure gradient may be neglected, 
and the velocity components are zero at the closed ends of the heat pipe. 

Also, his method of imposing a negligible radial pressure gradient by 
equating the mean axial pressure gradient over the cross section to the 
axial pressure gradient along the center line was used. 

Busse 7 began with the axial axisymmetric Navier-Stokes equation 
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where p is pressure, r is radial distance, u is radial velocity, 
v Is axial velocity, z is axial distance, n Is dynamic viscosity, and 
p is density. Busse Integrated the mass conservation equation with 
respect to u so that this term can be eliminated. This is also done in 
the present analysis except that density is variable in z and r: 
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When equation (2) is combined with equation (1) the following results: 
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Here, the axial variation of u has been neglected in the viscous term. 
The assumption is made that the radial variation of density is much less 
than the radial variation of the axial velocity since the latter vanishes 
at the wall. On the basis of this assumption, a presently undefined mean 
density p m and its axial derivative are extracted from the integrals, 
yielding 
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The crux of Busse's attack on the Incompressible case is the assumption 
of an approximate velocity profile which satisfies all the boundary con- 
ditions. It incorporates a parameter A(z) which may be constant or may 
vary with position, depending on mass influx boundary conditions. The 
same profile will be employed in this analysis: 


v = 




1 + A(z) 




(5) 
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Here, v m Is the mean axial velocity and R is the vapor space diam- 
eter. When A(z) = 0, equation (5) reduces to the Poiseuille profile. 
Integration of equation (4) over the cross section of the pipe and 
manipulation give 
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The mean density is assumed to be related to the pressure and a mean tem- 
perature T m by the ideal gas law, which yields upon differentiation 
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and where Rg, the gas constant, may depend on pressure and temperature 
if chemical reactions are considered. 


Continuity requires that the vapor mass flux m(z) at any cross section 
equals the product of mean density p m , the vapor cross-sectional area 
A v , and the mean axial velocity v m (z). From this relationship and 
equation (7) the axial derivative of v m is found to be: 
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Equations (6) to (8) are combined and the indicated integrations are per- 
formed. Dropping the distinction between the mean cross section pressure 
and the local pressure by Bus se 1 s approximation gives 
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where 
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The variation of A v with z is understood to be small. 

The mean temperature T m is evaluated by means of an energy integral. 
Figure 1 shows the control volume through which energy is transported. 

The energy flows into the boundary are equated with the energy flews out, 
internal energy and pressure are eliminated by using enthalpy, and the 
limit Az -» 0 i s taken, to give 
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where h is the enthalpy of the mixture per unit mass, qR is the energy 
flux entering from the wall of the pipe, and <l> is the viscous dissipa- 
tion per unit volume. The relation for ur, the wall radial velocity, is 
given by the mass conservation equation in the form 
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The energy influx q^ is specified as -p^u^Ce^ + u|/2) , where eR is 
internal energy. The result of combining equations (10) and (11) Is 
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The mean value of both the enthalpy and the density can be separately 
extracted from the integrals if the assumption is again made that the 
variation of both variables with radius is much less than that of axial 
velocity. When these mean values have been extracted, and terms common 
to all expressions cancelled, the energy equation becomes 
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If chemical equilibrium In the vapor must be considered, enthaloy h m 
becomes a function of both T m and p. When its derivative is inserted 
into equation (13) with equations (7) and (9), and the integrations are 
performed, the result is 
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An additional equation required to find the three dependent variables p, 
A, and T m comes from Busse's? approximation that the axial gradient of 
the mean pressure equals the center line pressure gradient. The latter, 
from equation (1) for r = 0 using equations (4) and (5), is 
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The mass flux is related to the heat flux into the pipe by 

dq R 1 dm [ h ( 1 dm\ 2 

dz " 2-rrR dz v/ + V2 ttR Pr dz J 


(16) 


where the kinetic energy of the injected vapor has been incorporated, 
and h V( r is the heat of evaporation. Equations (9), (14), and (15), 
with equation (16) and appropriate boundary conditions enable a solution 
for p, A, and T m . 
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The equations derived here pertain to alkali metal working fluids if 
equilibrium between monomeric and dimeric species is experienced. With 
proper application, they also apply to working fluids of a single gaseous 
species. If, in the case of the alkali metals, composition remains fixed 
from the point of injection, the equations must be modified, ihis case 
has been considered but the results are much less satisfactory than r or 
equilibrium and are not presented here. 


RESULTS AND DISCUSSION 

To illustrate the application of the algorithm, several cases were run 
Involving alkali liquid metal working fluids. Entropies and enthalpies 
of the species referred to a common base that includes chemical energy 
were used with polynomials for the specific heats to find thermodynamic 
properties. 3 These data enabled the determination of vapor pressure, 
heat of vaporization, various first derivatives and composition of the 
vapor phase. The latter includes both monomers (atoms) and dimers (mole- 
cules) in the case of alkali metals. The assumption is made that the 
ideal gas law applies to each of the vapor species. WhUe it would be 
desirable to incorporate real gas effects in the form of viria! equa- 
tions, this would greatly complicate the computations. 

The three equations involving the first derivatives of p,_A, and T m 
were solved by 3. simple Runge-Kuttei routine^ ror the set of s 1 mu 1 tdneous 
equations (9), (14), and (15). No study of the effect of step size on 
the precision of results was made. Each computation began with 20 uni- 
form steps in the evaporator, but step size was reduced as the slope of 
p became Increasingly steep in the vicinity of the sonic limit. 

Figure 2 illustrates a simple application of the algorithm which solves 
the equation set. The heat input to a sodium-filled heat pipe of 0.02 m 
diameter with an evaporator length of 0.1 m was increased until the pine 
was choked at an evaporator endcap vapor temperature of 900 K and a_heat 
input of 5870 W. The sonic limit computed from the approximation of 
Busse 10 was about 5580 W. Figure 2 shows the profiles of vapor mean tem- 
perature T m , wall temperature Tr, the velocity profile parameter A, 
and the pressure ratio p/po through the evaporator, where subscript 
0 denotes the vapor condition at z = 0, the endcap. Figure 3 shows the 
velocity profile v/v m at the downstream end of the evaporator where 
choking occurs, with a Poiseuille profile for comparison. The maximum 
value did not occur at the centerline, as expected. Busse's apciroximate 
profile is a three-point fit, so the profile shown in Figure 3 must be 
regarded as the best representation of a very flat profile available 
under the circumstances. 

In the next example two evaporator sections in sequence were considered, 
the first 0.2 m long, the second 0.1 m long, separated by an adiabatic 
section of 0.2 m length. The working fluid was sodium and the evaporator 
end cap temperature was 1000 K. Both evaporator sections were receiving 
heat at the rate of 71 600 W/m, which resulted in choke at the end of the 
second evaporator. Figure 4 shows the profiles of temperature, pressure 
ratio, and A. The scale of the figure is too small to show clearly tne 
decline in the parameter A in the adiabatic section from 0.738 to 
0.704. In i ncompress i bl e flow, the limit for A in very long adiabatic 
sections should decline toward zero, the Poiseuille condition. 
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Figure 2. - Variation of heat pipe parameters with location. 
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Figure 3. - Comparison of V/V m for Poisseiulle and chocked flow conditions. 
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Figure 4. - Variation of heat pipe parameters with location for a multi-location heat source on 
evaporator. 
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Another example shows the sonic limit heat transport against evaporator 
endcap temperature tor a scdi um-f i 1 1 ed evaporator of 0.02 m diameter ana 
0.1 m length. The short evaporator length minimized the viscous pressure 
loss, so the results can be comoared to those computed from inviscid 
equations of Busse 10 and Levy. 1,1 A simple control structure appended to 
the vapor pressure subroutine generated figure 5 in a short time. 

The most suitable data against which to check the algorithm are tnose of 
Ivanovskii et al. 6 They inserted a temperature microprobe into a sodium 
heat pipe having an evaporator length of 0.1 m, an adiabatic section 

length of 0.05 m, a condenser length of 0.55 m, and a diameter of 

0.014 m. Temperatures were read at intervals along the axis with the 
probe. Liquid-saturated wicking on the probe was to ensure that the tem- 
peratures being observed were the same as those at the vapor-liquid 
interface. No details were given concerning any effect of the pressure 
field around the probe tip on the measured temperature. For an input of 
about 500 W, stated to be accurate to within 6 to 10 percent, tne data 

are shown in figure 6. Also shown are lines of T^ and Tq computed 

by the algorithm. Equations (9), (14), and (15) were used only in the 
evaporator and adiabatic sections. An extension into the condenser was 
made by equating T m to Tq beyond the point where T m exceeded Tq. 
The justification is the observation that at seme point in the concenser, 
Trp must converge to Tq as it did at the beginning of the evaporator. 
For this portion of the solution, the equations were suitably modified 
to yield as dependent variables only p and A. This simple procedure, 
enabling the solution to be continued into the condenser, neglected the 
known complications of condensation shocks and recirculation. 

For the same heat pipe at a heat input of 1000 W, figure 7 shows experi- 
mental data of Ivanovskii and lines of T m and Tq from the algorithm. 
Also shown is a curve of Tq computed by Chen and Faghri 1 ^ using tne 
Phoenics code. It is in basic agreement with the results computed by 
the algorithm presented here. The Phoenics code, a sophisticated finite 
difference procedure, can handle far more complex situations than the 
elementary algorithm herein. These include evaporators with very small 
ratios of length to diameter and heat pipes of complex cross section. 

The size of computer required for the Phoenics code and its cost make it 
less suitable for everyday heat pipe analysis. Furthermore, the simple 
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TEMPERATURE, K 
Figure 5. - Sonic limit determination. 
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Figure 6. - Comparison of predicted values with ivanovskii data. 



Figure 7. - Comparison of predicted values with Ivanovskii data and results of other methods. 

equations derived from first principles in the present caper may facili- 
tate a better physical understanding of the governing phenomena. 


CONCLUDING REMARKS 

The heat pipe vapor flow algorithm discussed herein has been shewn to 
model the compressible vapor dynamics in the evaporator and aaiabatic 
sections of pipes of conventional proportions. The model has introduced 
complexities which have usually been ignored. These include chemical 
equilibrium in the case of liquid metals, a velocity profile parameter 
dependent upon pressure and temperature, and a mean vapor temperature 
differing from the vapor interface temperature. Future work might be 
expended profitably on the introduction of a reasonaole temperature pro- 
file into the energy equation, and the inclusion of droplet nucleation 
rates into the momentum and energy equations. Also, a satisfactory 
representation of the vapor effects in the condenser requires that the 
transition to turbulence, the treatment of condensation shocks, flow 
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recirculation , and other phenomena be modeled concisely for use in simple 
programs on small computers. 
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